/** Fichier implémentant un objet donné par un potentiel en @f$ r^2 @f$.
 *
 * Ce fichier contient l'implémentation de la classe MetaBalle.
 *
 * @file MetaBalle.cpp
 * @author Johan "Solo" Luisier
 * @date 27 Février 2013
 */
#include "MetaBalle.hpp"

#include "Algebre/ConstruitMatrice.hpp"
#include "Algebre/Solveur.hpp"

#include "Utils/ATL.hpp"

using namespace std;

MetaBalle::MetaBalle( const vector< Coord >& poles,
                      const Vecteur& xPrime,
                      const Vecteur& yPrime,
                      const double& portee,
                      const double& echelle )
    : Equipotentielle( "MetaBalle" )
    , LongueurInteraction( portee )
    , FacteurEchelle( echelle )
    , Poles( poles.begin(), poles.end() )
{
    MatriceSO3 mat;
    StatusCode sc( axesXY( xPrime, yPrime, mat ) );

    if ( ! sc )
    {
        cerr << sc.message() << endl;
        exit( 1 );
    }

    if ( poles.empty() )
    {
        cerr << "Aucun pôle donné pour MetaBalle" << endl;
        exit( 1 );
    }

    vector< Coord >::const_iterator polIt, polFin( poles.end() );

    for ( polIt = poles.begin(); polIt != polFin; polIt++ )
        PseudoCentre += *polIt;

    PseudoCentre /= double( poles.size() );

    VersLocal = Transformation( mat, - PseudoCentre, 1. / FacteurEchelle );
    VersGlobal = VersLocal.inverse();

    PseudoRayon = 2. * LongueurInteraction + 4.;

    MoteurResolution = new Solveur();
}

MetaBalle::~MetaBalle()
{}

Coord MetaBalle::trouveIntersection( const Rayon &ray ) const
{
    Coord pt( VersLocal( ray.origine() ) );
    Vecteur dir( VersLocal( ray.direction() ) );

    /**
     * On calcule l'intersection avec l'équation pour @f$ r < \frac{L}{3} @f$.
     */
    double a( 0. ), b( 0. ),
            c( ( PotentielDeSurface - /*FacteurEchelle **/ double( Poles.size() ) ) * puissance< 2u >( LongueurInteraction ) );

    vector< Coord >::const_iterator polIt, polEnd( Poles.end() );

    for ( polIt = Poles.begin(); polIt != polEnd; polIt++ )
    {
        a += dir.norme();
        b += 2. * dir * ( pt - VersLocal( *polIt ) );
        c += ( pt - VersLocal( *polIt ) ).r2();
    }

    MoteurResolution -> setParams( c, b, a );

    vector< double > solutions;

    MoteurResolution -> resoud( solutions );

    vector< double >::const_iterator solIt, solFin( solutions.end() );

    Coord local( false );

    for ( solIt = solutions.begin(); solIt != solFin; solIt++ )
    {
        if ( *solIt > 1.e-7 )
        {
            local = pt + (*solIt) * dir;
            break;
        }
    }

    //cerr << "local = " << local << ", norme = " << local.r() << endl;

    return local;
}

VecteurU MetaBalle::normale(Coord const &pt) const
{
    Coord local( VersLocal( pt ) );

    double dist( local.r() );

    Vecteur loc;

    if ( dist < LongueurInteraction )
        loc = 2. /* * FacteurEchelle*/ / puissance< 2u >( LongueurInteraction ) *
                Vecteur( pt.x() , pt.y(), pt.z() );

    return VecteurU( VersGlobal( loc ) );
}

/** Permet de vérifier si un point est à @e l'intérieur de l'instance.
 *
 * On vérifie que le potentiel du point testé est @e supérieur au potentiel de
 * la surface.
 *
 * @param[in] pt point à tester
 *
 * @retval true si le point est bien dans le volume, @retval false sinon.
 */
bool MetaBalle::estDans(const Coord & pt) const
{
    return potentiel( pt ) - PotentielDeSurface > -1.e-7;
}

/** Permet de vérifier si un point est à la @e surface de l'instance.
 *
 * On vérifie que le potentiel du point testé est @e égal au potentiel de
 * surface.
 *
 * @param[in] pt point à tester
 *
 * @retval true si le point est bien sur la surface, @retval false sinon.
 */
bool MetaBalle::estSur(const Coord & pt) const
{
    return abs( potentiel( pt ) - PotentielDeSurface ) < 1.e-7;
}

void MetaBalle::mapping(const Coord &, int &i, int &j, const int &frequence, const int &Hauteur, const int &Largeur) const
{
    i = j = 0;
}


ostream& MetaBalle::affiche( ostream& os ) const
{
    return os << Nom << endl
              << "Volume centré en " << PseudoCentre << endl
              << "Pseudo rayon " << PseudoRayon << endl;
}

Coord MetaBalle::intrinseque( const Coord& ) const
{
    return Coord();
}

const double& MetaBalle::longueurInteraction() const
{
    return LongueurInteraction;
}

const double& MetaBalle::facteurEchelle() const
{
    return FacteurEchelle;
}

double MetaBalle::potentiel( const Coord& pt ) const
{
    /**
     * On passe en coordonnées locales pour évaluer le potentiel.
     */
    Coord local( VersLocal( pt ) );

    double resultat( 0. );

    double dist( local.r() );

    vector< Coord >::const_iterator polIt, polFin( Poles.end() );

    for ( polIt = Poles.begin(); polIt != polFin; polIt++ )
    {
        dist = ( local - *polIt ).r();

        if ( dist < LongueurInteraction )
            resultat += FacteurEchelle * ( 1. - puissance< 2u >( dist / LongueurInteraction ) );;
    }

    return resultat;
}
